Aprendizado Supervisionado

O objetivo desse notebook é aplicar a Linguagem R em Ciência de Dados na resolução de problemas de classificação.

Classificando flores

A base de dados iris fornece medidas em centímetros das variáveis largura e comprimento da sépala e largura e comprimento da pétala, respectivamente, para 50 flores, de cada uma das 3 espécies de íris.

data('iris')
head(iris)

Podemos usar o scatterplot para entender o comportamento das classes

plot(iris[,1:4], col=iris$Species)

É muito importante tratar valores ausêntes. A função is.na e is.nan nos diz se um elemento é NA ou NaN. Combinada com a função any, podemos descobrir se algum elemento da base é NA.

print(any(is.na(iris)))

## [1] FALSE

print(any(apply(iris, 2, is.nan)))

## [1] FALSE

Devemos normalizar os dados, uma vez que existem classificadores que são suscetíveis a escala dos atributos.

iris = data.frame(scale(iris[,1:4]), Species=iris$Species)
summary(iris)

##   Sepal.Length       Sepal.Width       Petal.Length      Petal.Width     
##  Min.   :-1.86378   Min.   :-2.4258   Min.   :-1.5623   Min.   :-1.4422  
##  1st Qu.:-0.89767   1st Qu.:-0.5904   1st Qu.:-1.2225   1st Qu.:-1.1799  
##  Median :-0.05233   Median :-0.1315   Median : 0.3354   Median : 0.1321  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.67225   3rd Qu.: 0.5567   3rd Qu.: 0.7602   3rd Qu.: 0.7880  
##  Max.   : 2.48370   Max.   : 3.0805   Max.   : 1.7799   Max.   : 1.7064  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

Depois de normalizar a base, precisamos particionar o conjunto de treinamento e teste. Para essa e outras bases vamos usar o holdout com 70% para treinamento e 30% para teste.

set.seed(123)
idx = sample(nrow(iris), 0.7*nrow(iris), replace=FALSE)
tran = iris[idx,]
test = iris[-idx,]

Podemos conferir se ocorreu contaminação.

intersect(rownames(tran), rownames(test))

## character(0)

Agora podemos induzir alguns modelos de Aprendizado de Máquina. Vamos começar pela Árvore de Decisão.

require('rpart')

## Loading required package: rpart

model = rpart(Species~., iris)

No caso da Árvore de Decisão, é possível plotar o modelo. Para isso vamos usar a biblioteca rattle.

library('rattle')

## Loading required package: tibble

## Loading required package: bitops

## Rattle: A free graphical interface for data science with R.
## Version 5.4.0 Copyright (c) 2006-2020 Togaware Pty Ltd.
## Type 'rattle()' to shake, rattle, and roll your data.

fancyRpartPlot(model, caption=NULL)

Também podemos predizer a classe das amostras de teste.

pred = predict(model, test, type='class')

E plotar o scatterplot para os dados de treinamento e teste.

par(mfrow=c(1,2))
plot(tran[,3:4], col=tran$Species, pch=as.numeric(tran$Species))
plot(test[,3:4], col=pred, pch=as.numeric(pred))

Podemos avaliar de forma mais objetiva o desempenho do modelo.

conf = table(test$Species, pred)
print(conf)

##             pred
##              setosa versicolor virginica
##   setosa         14          0         0
##   versicolor      0         18         0
##   virginica       0          1        12

Inclusive, calcular a taxa de acerto.

acc = sum(diag(conf))/sum(conf)
print(acc)

## [1] 0.9777778

Outras meedidas de desempenho poderiam ser utilizadas. Recomendo olhar os pacotes ‘caret’ e ‘mlr’ para mais informaões..

#F1-sconre
f1 <- function(conf) {
  prc = diag(conf)/colSums(conf)
  rec = diag(conf)/rowSums(conf)
  aux = (2 * prc * rec) / (prc + rec)
  return(mean(aux))
}

#AUC
require('ROCR')

## Loading required package: ROCR

auc <- function(class, pred) {
  aux = prediction(as.vector(pred), as.vector(class))
  aux = unlist(performance(aux, "auc")@y.values)
  return(aux)
}

Agora podemos induzir outros modelos de Aprendizado de Máquina. Para isso precisamos carregar algumas bibliotecas.

require('kknn')

## Loading required package: kknn

require('e1071')

## Loading required package: e1071

require('randomForest')

## Loading required package: randomForest

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:rattle':
## 
##     importance

Gerando os modelos de Aprendizado de Máquina.

model1 = kknn(Species~., tran, test, k=3)
model2 = svm(Species~., tran, kernel="radial")
model3 = randomForest(Species~., tran)

Calculando a predição.

pred1 = model1$fitted.values
pred2 = predict(model2, test, type='class')
pred3 = predict(model3, test)

Calculando o desempenho baseado na acurácia.

acuracia <- function(orig, pred) {
  acc = table(orig, pred)
  sum(diag(acc))/sum(acc)  
} 

cat('A acurácia do 3NN é', acuracia(test$Species, pred1))

## A acurácia do 3NN é 0.9333333

cat('A acurácia do SVM é', acuracia(test$Species, pred2))

## A acurácia do SVM é 0.9777778

cat('A acurácia do RF é', acuracia(test$Species, pred3))

## A acurácia do RF é 0.9777778

No entanto, o holdout não é a melhor alternativa. A validação cruzada estratificada é mais interessante. Podemos ver na imagem um exemplo de seu funcionamento.

k-fold cross validation (Adaptado da Wikepedia)

Podemos implementar :)

kfolds <- function(y, folds) {

  names(y) =  1:length(y)
  index = lapply(1:nlevels(y), function(i) {
    rep(1:folds, length.out=length(y[y == levels(y)[i]]))
  })

  index = unlist(index)
  folds = lapply(1:folds, function(i) {
    as.numeric(names(y[index == i]))
  })

  return(folds)
}

E agora podemos usar.

idx = kfolds(iris$Species, 10)

Antes vamos conferir se ocorreu contaminação.

intersect(idx[[1]], idx[[2]])

## numeric(0)

E agora, de fato, podemos usar.

acc = lapply(idx, function(i){
  model = randomForest(Species~., iris[-i,])
  pred = predict(model3, iris[i,])
  acuracia(iris[i,]$Species, pred)
})

cat(mean(unlist(acc)), '+-', sd(unlist(acc)))

## 0.9933333 +- 0.02108185

Classificando dígitos

Um experimento interessante é a tarefa de classificar dígitos utilizando o algoritmo de Redes Neurais Artificiais. Para isso precisamos relizar os experimentos em duas fases: (1) treinar a rede para reconhecer os dígitos e (2) utilizá-la para classificar novos dígitos. Para isso precisamos carregar as blbiotecas ‘tensorflow’ e ‘keras’.

require('tensorflow')

## Loading required package: tensorflow

require('keras')

## Loading required package: keras

A base de dados de dígitos se chama MNIST e pode ser carregada com a função dataset_mnist. Ela contém 60000 amostras dos dígitos 0 até 9. Cada dígito é uma imagem de 28 x 28 pixels.

mnist = dataset_mnist()

A função show_digit recebe um dígito na forma de uma matriz 28 x 28 e imprime o dígito na forma de imagem.

show_digit <- function(numero) {
  image(numero[1:28,28:1])
}

Podemos utilizar a função show_digit para plotar uma amostra de todos os dígitos presentes na base de dados.

par(mfrow=c(2,5))
plot = lapply(0:9, function(i) {
  aux = mnist$train$x[mnist$train$y == i,,][1,,]
  show_digit(t(aux))
})

Normalizando os dados de treinamento e teste.

mnist$train$x = mnist$train$x/255
mnist$test$x = mnist$test$x/255

Construindo uma Rede Neural Artificial usando Keras. Essa rede é composta de uma camanda de entrada flatten (achatada de 28 x 28 neurônios), uma camada oculta dense (128 neurônios com função de ativação linear) e uma camada de saída dense (10 neurônios com função de ativação exponencial normalizada).

modelo = keras_model_sequential() %>% 
  layer_flatten(input_shape = c(28, 28)) %>% 
  layer_dense(units = 128, activation = "relu") %>% 
  layer_dense(10, activation = "softmax")

Informações sobre a Rede Neural Artificial.

summary(modelo)

## Model: "sequential"
## ________________________________________________________________________________
## Layer (type)                        Output Shape                    Param #     
## ================================================================================
## flatten (Flatten)                   (None, 784)                     0           
## ________________________________________________________________________________
## dense_1 (Dense)                     (None, 128)                     100480      
## ________________________________________________________________________________
## dense (Dense)                       (None, 10)                      1290        
## ================================================================================
## Total params: 101,770
## Trainable params: 101,770
## Non-trainable params: 0
## ________________________________________________________________________________

Podemos representar graficamente:

library(deepviz)
library(magrittr)
modelo %>% plot_model()

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Ainda precisamos definir o algoritmo de treinamento, a função de avaliação e a medida de desempenho que vamos utilizar. Podemos fazer isso com a função compile.

modelo %>% 
  compile(
    loss = "sparse_categorical_crossentropy",
    optimizer = "adam",
    metrics = "accuracy"
  )

Agora podemos treinar a rede.

historico = modelo %>% 
  fit(
    x = mnist$train$x, y = mnist$train$y,
    epochs = 10,
    validation_split = 0.3,
    verbose = 2
  )

É interessante medir o erro e o desempenho da rede ao longo das épocas no conjunto de treinamento e validação.

plot(historico)

## `geom_smooth()` using formula 'y ~ x'

Agora que sabemos que a rede aprendeu a reconhecer dígitos, podemos verificar seu desempenho no conjunto de teste. Para isso precisamos aplicar essas amostras nunca antes vistas na entrada da rede e obter os valores de saída.

predicao = predict(modelo, mnist$test$x)
head(predicao)

##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,] 1.426630e-09 5.901615e-11 3.042254e-06 7.273457e-05 3.330658e-12
## [2,] 2.037992e-09 2.025005e-08 9.999962e-01 3.677915e-06 2.707853e-20
## [3,] 6.741865e-07 9.969687e-01 5.369112e-04 1.404337e-05 1.767417e-06
## [4,] 9.999970e-01 2.808441e-12 3.241617e-07 7.545845e-11 1.836492e-08
## [5,] 1.148857e-09 9.923532e-13 1.644524e-08 6.411225e-12 9.999746e-01
## [6,] 6.497803e-10 9.994200e-01 2.336848e-06 5.812358e-08 2.098802e-08
##              [,6]         [,7]         [,8]         [,9]        [,10]
## [1,] 1.024663e-08 1.086371e-15 9.999231e-01 3.681009e-09 1.066019e-06
## [2,] 3.542250e-09 4.432574e-09 1.772909e-13 1.003974e-07 1.799808e-17
## [3,] 2.197488e-06 2.226117e-05 1.876961e-03 5.764504e-04 2.067261e-07
## [4,] 3.680925e-08 9.648173e-07 2.516348e-07 4.885596e-11 1.386067e-06
## [5,] 4.852368e-12 1.236614e-08 4.410211e-07 9.589456e-10 2.486780e-05
## [6,] 8.363623e-11 3.038747e-09 5.757674e-04 1.727228e-06 3.193629e-10

Cada neurônio de saída da rede irá retornar um valor. Aquele neurônio com maior valor deve indicar a classe com maior proababilidade.

predicao = apply(predicao, 1, which.max)
predicao[1:6]

## [1] 8 3 2 1 5 2

Calculando o desempenho da rede no conjunto de teste.

acc = acuracia(mnist$test$y, predicao)
cat('A acurácia do modelo no conjunto de teste eh:', acc)

## A acurácia do modelo no conjunto de teste eh: 0.9733

Agora podemos usar os modelos de Árvore de Decisão, kNN, SVM e Random Forest para comparar com o desempenho com a Rede Neural Artificial. Primeiro vamos construir uma base tabular com a mnist.

tran = data.frame(mnist$train$x, y=as.factor(mnist$train$y))
test = data.frame(mnist$test$x, y=as.factor(mnist$test$y))

Agora podemos induzir outros modelos.

model = rpart(y~., tran)
pred = predict(model, test, type='class')

cat('A acurácia da ANN é:', acc)

## A acurácia da ANN é: 0.9733

cat('A acurácia da DT é', acuracia(test$y, pred))

## A acurácia da DT é 0.6196